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Abstract 

Testing the constancy of the gravitational constant G is a longstanding fundamental question in 
natural science. As first suggested by Jofre, Reisenegger and Fernandez [Ij], Dirac's hypothesis of a 
decreasing gravitational constant G with time due to the expansion of the Universe would induce 
changes in the composition of neutron stars, causing dissipation and internal heating. Eventually, 
neutron stars reach their quasi-stationary states where cooling, due to neutrino and photon emis- 
sions, balances the internal heating. The correlation of surface temperatures and radii of some old 
neutron stars may thus carry useful information about the rate of change of G. Using the density 
dependence of the nuclear symmetry energy, constrained by recent terrestrial laboratory data on 
isospin diffusion in heavy-ion reactions at intermediate energies, and the size of neutron skin in 
within the gravitochemical heating formalism developed by Jofre et al. ]|, we obtain an 
upper limit for the relative time variation \G/G\ in the range (4.5 — 21) x 10^^'^yr~^ . 
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I. INTRODUCTION 



The question whether or not the fundamental constants of nature vary with time has been 
of considerable interest in physics. The constancy of the gravitational coupling parameter 
G was first addressed in 1937 by Dirac [2] who suggested that the gravitational force might 
be weakening due to the expansion of the universe. Although general relativity assumes a 
strictly constant G, time variations of the Newton's constant are predicted by_^some alter- 

|6[[. Many 



native theories of gravity fsl and a number of modern cosmo 
theoretical approaches, such as models with extra dimensions 

, Q, Q, H, Q 



and scalar-tensor quintessence models 



11 



12 



7|, string theories [8| 



I, a 



id, 



have been proposed in which 



the gravitational coupling parameter becomes a time- dependent quantity. Nowadays the de- 
bate over the constancy of G has been revived by recent astronomical observations [2^ 
of distant high-red-shift type la supernovae suggesting that presently the Universe is in a 
state of accelerated expansion [6]. This acceleration can be interpreted in terms of a "dark 
energy" with negative pressure, or alternatively by allowing a time variation of the gravita- 



tional constant 
and Kothari 



J]. Soon after Dirac had published his hypothesis |2|], Chandrasekhar 17 1 



181 ] pointed out that a decreasing G with time could have some detectable as- 



trophysical consequences. Since then many attempts have been made to find astrophysical 
signs due to the possible time variation of G. However, there is no firm conclusion so far (see 
Ref. [2IJ for a review). Interestingly though, as pointed out by Uzan 2l|], contrary to most 
of the other fundamental constants, as the precision of the measurements increased the dis- 
crepancy among the measured values of G also increased. This circumstance led GOD AT/. 
(Committee on Data for Science and Technology) to raise the relative uncertainty for G 21] 
by a factor of about 12 in 1998. Some of the previous upper limits on the time variation 
of G, as obtained by different experiments/methods, are summarized in Table 1 (adapted 
from Reisenegger et al. [69i])- Given the current status of both theory and experiment, it is 
fair to say that whether or not the gravitational constant varies with time is still an open 
question and therefore additional work is necessary to investigate further this fundamental 
issue. 

Recently a new method, called gravitochemical heating jl|, has been introduced to con- 
strain a hypothetical time variation in G, most frequently expressed as \G/G\. In Ref. jl| the 
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TABLE 1: Upper bounds on the time variations of G (adapted from Reisenegger et al. [691] ). 
The first column hsts the method used. The second column contains the upper limit on the time 
variation of G, most usefully expressed as |G/G|, the third column is a rough time scale over which 
each experiment is averaging this variation, and the last column is the corresponding reference. 
The first two experiments on the list probe the variation of G from the early Universe to the present 
time. The next four experiments are sensitive to long time-scales, but without reaching to the early 
Universe. And the last three experiments on the list probe the change of G over short time-scales 
of years and decades. 



Method 


|G/GU,,,[10-i2y^-i] 


Time scale [yr] Reference 


Big Bang Nucleosynthesis 


0.4 


1.4 X 10^° 


[70] 


Microwave Background 


0.7 


1.4 X 10^° 


[7JJ 


Gloubular Cluster Isochrones 


35 


IQio 


[72] 


Binary Neutron Star Masses 


2.6 


IQio 


m 


Helioseismology 


1.6 


4 X 10^ 


m 


Paleontology 


20 


4 X 10^ 


[75] 


Lunar Laser Ranging 


1.3 


24 


[7^] 


Binary Pulsar Orbits 


9 


8 


m 


White Dwarf Oscillations 


250 


25 


[78] 



authors suggested that such a variation of the gravitational constant would perturb the in- 
ternal composition of a neutron star, producing entropy which is partially released through 
neutrino emission, while a similar fraction is eventually radiated as thermal photons. A 
constraint on the time variation of G is achieved via a comparison of the predicted surface 
temperature with the available empirical value of an old neutron star [22]. The gravito- 



2J] (see 



chemical heating formalism is based on the results of Fernandez and Reisenegger 
also who demonstrated that internal heating could result from spin-down compression 
in a rotating neutron star {rotochemical heating). In both cases (gravito- and rotochemical 
heatings) predictions rely heavily on the equation of state (EOS) of stellar matter used to 
calculate the neutron star structure. Accordingly, detailed knowledge of the EOS is critical 
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for setting a reliable constraint on the time variation of G. The global properties of neutron 
stars such as masses, radii, moments of inertia, thermal evolution, etc have been studied 
extensively, see, e.g.. 



26 



27 



28 



29 



30 



31 



32 



34| . Generally, predictions differ widely. 



mainly, due to the uncertainties of the equations of state employed in neutron star structure 
calculations [34]. Therefore, determining the EOS of stellar matter is a question of central 
importance with an answer requiring understanding of proper nuclear physics. 

Using the well tested parabolic approximation the EOS of isospin asymmetric nuclear 
matter is written as 

e(p, a) = e(p, 0) + esym(p)a^ (1) 

which manifests the separation of the EOS into isospin symmetric (the energy per particle 
of symmetric nuclear matter) and isospin asymmetric contributions (the nuclear symmetry 
energy egym)- In the above expression a = {pn — Pp)/p is the usual asymmetry parameter, 
P = Pn + Pp the baryon number density, and pn and Pp are the neutron and proton densities 
respectively. While the EOS of symmetric nuclear matter {a = 0) is relatively well under- 
stood, the density dependence of the symmetry energy, egym, is still very poorly constrained 
especially at high densities. Variations in the esym predicted by various models often yield 
dramatically different predictions for properties of neutron stars (e.g., see Ref. 3J]). Be- 
cause of its importance for neutron star structure, determining the density dependence of 
the symmetry energy has been a high-priority goal for the intermediate energy heavy-ion 
community. Although extracting the symmetry energy is not an easy task due to the compli- 



cated role of isospin in reaction 



have been suggested 



35 



36 



37 



dynamics, several promisin g pr obes of the symmetry energy 



38| (see also Refs. [39|, |40|, |4l| for reviews) 



Some significant progress has been made recently in determining the density dependence 
of Csym using: (1) isospin diffusion in heavy-ion reactions at intermediate energies as a probe 
of both the magnitude and slope of the symmetry energy around the saturation density 

(2) flow in heavy-ion collisions at higher energies to 



42 



43 



44 



451, 



46 



(po ~ 0.16/m 

constrain the equation of state of symmetric nuclear matter 401], and (3) the sizes of neutron 



skins in heavy nuclei to constrain egymip) at subsaturation densities 



32 



48|,|49| 



In this work, we combine recently obtained diffusion data, information from flow observ- 



ables, studies on the neutron skin of "^^^Pb, and other information to constrain a possi 
time-variation of the gravitational constant G through the gravitochemical formalism 
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We do not aim to add anything fundamental to the original method of Ref. [l|]. Our ob- 
jective is to provide a restrictive upper limit for the time variation of G, applying an EOS 
constrained by terrestrial empirical data from nuclear reactions induced by neutron-rich nu- 
clei. After the introductory notes in this section, we discuss, in some details, the general 
formalism of the gravitochemical method. The equations of state used in this study are 
outlined briefly in Section III. Our results for the upper limits of \G/G\ are presented and 
discussed in Section IV. The effects of "exotic" (hyperonic/quark) phases in neutron star 
matter on a possible time variation of G are addressed in Section V. We conclude in Section 
VI with a short summary. 



II. GRAVITOCHEMICAL HEATING: FORMALISM 

To provide the reader with a self-contained manuscript, in this section we recall the 
main steps leading to the calculation of the surface temperature of an old neutron star via 

nn 

the gravitochemical heating method. For a detailed discussion see Refs. [11,124]. (Conven- 
tions and notation as in the above references.) The simplest neutron star models assume a 
composition of nucleons and light leptons, electrons and muons, which can transform into 
each other through direct and inverse /3-reactions. The neutrinos {v) and antineutrinos {v) 
produced in these reactions leave the star without further interactions, contributing to its 
cooling. In /3-equilibrium the balance between the rates of direct and inverse processes is 
reflected through the following relation among the chemical potentials of the particle species 

Hn- l^p = fJ'e = fJ'i, (2) 

I 1 

As pointed out in Ref. [!)] a time-variation of G would cause continuously a perturbation 
in the stellar density and since the chemical potentials are density-dependent the star thus 
always departs from /3-equilibrium. This departure is quantified by the chemical imbalances 

Vnpe = Sfin " Sfip - 6fle (3) 
Vnpfj, = Sfin - Sfip - djjf, (4) 

where SfXi = /ij — /i^^ is the deviation of the chemical potential of particle species i {i = 
n, p, e, /i) from its equilibrium value at a given pressure. The chemical imbalances enhance 
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the rates of reactions driving the star to a new equihbrium state. If G changes continuously 
with time the star will be always out of equilibrium, storing an excess of energy that is 
dissipated as internal heating and enhanced neutrino emission 

The evolution of the internal temperature is given by the thermal balance equation 

T^ = ^[Lh-LT-L^] (5) 

where C is the total heat capacity of the star, LJ^ is the total power released by heating 
mechanisms, is the total neutrino luminosity, and L'^ is the photon luminosity ("oo" 
labels the quantities as measured by a distant observer). The evolution of the red-shifted 
chemical imbalances is governed by 

= ^^^n - Sfl^ - SfiT (6) 
C = Sf^n - Sf^p - (7) 

These equations can be written as |l| 



Inpe 

oo 

'npfj, 



= - [AnAVupe^ + AMAVnpe. T^] " [BdAVup,^ T^) + BmAVup,^ T^)] (8) 
= - [AnAVnpe, Tn + AmAvTps, T^] - [BdAVup,, T^) + BmAvTp,, T^)] (9) 

The functions A and B quantify the effect of reactions toward restoring chemical equilibrium, 
and thus have the same sign as rjnpi (/ = e,/i) [24]. The subscripts D refers to the so-called 
direct Urea cooling processes 

n — >■ p + I + u 
p + l n + u, (10) 

which are fast but possibly forbidden by the energy-momentum conservations when the 
proton fraction is low The subscripts M refers to the modified Urea reactions 

n + N ^ p + N + 1 + 9 
p + l + N n + N + v, (11) 

which are slow and an additional nucleon or nucleus must participate in order to con- 
serve momentum |29|. The constants Cnpe and Cnpfi quantify the departure from chemical 
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equilibrium due to a time-variation of G. They can be written as jl| 

Cnp^ — {Znpn Zfip) Iq^^ -\- ZfiplQ p (12) 

Here I^i = {ON^'^ / dG)A is the change of the equihbrium number of particle species i {i = 
n,p, e, /i), iVj^''', due to the variation of G and Z are constants depending only on the stellar 
structure |l!]. Equations (5), (8), and (9) determine completely the thermal evolution of 
a neutron star with gravitochemical heating. The main consequence of this mechanism 
is that eventually the_ star arrives at a quasi-equilibrium state, with heating and cooling 
balancing each other [1]. The properties of this stationary state can be obtained by solving 
simultaneously Eqs. (5), (8), and (9) by setting T°° = rj'^^ = t)^^ = 0. The existence 
of a quasi-equilibrium state makes it possible, for a given value of \G/G\, to compute the 
temperature of an old neutron star without knowing its exact age [1], since, due to the 
independence of the solution from the initial conditions, it is unnecessary to model the 
complete evolution of the chemical imbalances and temperature. 

If only the modified Urea reactions are allowed, for a given stellar model, it is possible to 
derive an analytic expression relating the photon luminosity in the stationary state, L^^g, to 
\G/G\. This is because the longer time scale required to reach a stationary state, when onl 
modified Urea processes operate, results in chemical imbalances satisfying rjnpi >> ksT 
Under these conditions the photon luminosity in the quasi-equilibrium state is given by 



nly 



LT.. = Cm 



c 



H 



t8 \ 1/7 / r8 \ 1/7- 



8/7 

(13) 



The meaning of the constants Cm and Ch, and the functions Lm^ {i = e,/i) is explained in 
Refs. . From L^^eq ^^e neutron-star surface temperature can be calculated by assuming 

an isotropic blackbody spectrum 

L-^ = 47rcri?^(rr)^ (14) 

with a the Stefan-Boltzmann constant and i?oo the red-shifted radius of the star. In the 
case when only slow /5-reactions operate, we write the stationary surface temperature as 

A 2/7 

^ G ' (^^^ 
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where the function P is a quantity depending only on the stellar model and the equation 
of state. On the other hand, if at higher densities the proton fraction is large enough so 
that the direct Urea reactions are allowed, the thermal evolution of a neutron star with 
gravitochemical heating needs to be modeled by solving numerically the coupled Eqs. (5), 
(8) and (9). 

As demonstrated by Jofre et al. [ll the formalism outlined here can be applied to constrain 
the value of \G/G\, provided one knows (i) the surface temperature of a neutron star, and 
(ii) that the star is certainly older than the time-scale necessary to reach a quasi-stationary 
state. So far the only object satisfying both conditions is PSR J0437-4715, which is the 
closest millisecond pulsar to our solar system. Its surface temperature was deduced from 
ultraviolet observations 22] while its mass was determined by Hotan et al. [3] to be in 
the range MpsR = (1.1 — 1.5)Mq. (Another mass constraint, MpsR = 1.58 ± 0.18Mq, was 
given previously by van Straten et al. \2^.) To constrain the value of \G/G\ one, therefore, 
needs to consider neutron-star models in the above mass range and calculate the surface 
temperature for each stellar configuration. 



III. EQUATION OF STATE AND NEUTRON STAR STRUCTURE 

Clearly, predictions of the surface temperature and, in turn, value of \G/G\ depend heavily 
on the EOS of neutron-star matter since the later is critical for determining the neutron- 
star structure. Currently, theoretical predictions of the EOS of neutron-rich matter diverge 
widely mainly due to the uncertain density dependence of the nuclear symmetry energy. 
Consequently, to provide a stringent constraint on the time variation of G, one should 
attempt to reduce the uncertainty due to the esym{p)- Recently available nuclear reaction 
data allowed us to constrain significantly the density dependence of the symmetry energy 
mostly in the sub-saturation density region. While high energy radioactive beam facilities 
under construction will provide a great opportunity to pin down the high density behavior 
of the nuclear symmetry energy in the future. In this work, we apply the gravitochemical 
method with several EOSs describing matter of purely nucleonic (npe/i) as wells as hyperonic 
and hybrid stars. Among the nucleonic matter EOSs, we pay special attention to the one 

energy esym{p) of the MDI EOS 
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is constrained in the sub-saturation density region by the available nuclear laboratory data, 
while in the high-density region we assume a continuous density functional. The EOS of 
symmetric matter e(p, 0) for the MDI interaction is constrained up to about five times the 
normal nuclear matter density by the available data on collective flow in relativistic heavy- 
ion reactions. 

Let us first briefly recall here the main ingredients of the MDI EOS following Ref. 52 1. 
The MDI EOS corresponds to the single-particle potential 



U{p,a,p,T,x) = A^{x)^ + Mx)^ + B[^] (1 



Pr 



P 



xa 



Po 



+ 



Po 



.Po, 



+ 



3 / 



dy 



B p-i 

QTX ap-r 

^+1 Po 



(16) 



deduced 



53l | from the Gogny interaction. In the above equation x is a parameter introduced 



to reflect the largely uncertain density dependence of the Csymip) as predicted by various 
many-body approaches; t{t') is 1/2 (—1/2) for neutrons (protons) with t ^ t'\ o = 4/3, 
fri^yP) is the space distribution function at coordinate r and momentum p; A^, Ai, B, Cr^r, 
Cr'^r' and A are parameters fixed by fitting the momentum dependence of U{p, a,p, r, x), as 
predicted by the Gogny/Hartree-Fock and/or Brueckner-Hartree-Fock (BHF) calculations, 
so that the saturation properties of nuclear matter and the value of the symmetry energy 
at the saturation density (esym(po) ~ 32MeV) are predicted correctly. The compression 
modulus of saturated nuclear matter, k, is set to 211MeV consistent withe empirical range 



recently suggested by Garg 



5|]. More specifically, B = 106.35MeV and A = A;^ is the 



nucleon Fermi momentum in symmetric nuclear matter. The quantities A„(x) and Ai{x) 
depend on the parameter x according to 

2B 



AJx) 



Aiix) 



-95.98 -X- 



-120.57 



cr + 1 
2B 

X- 



a 



1 



(17) 
(18) 



The isoscalar potential, evaluated from Uq = {Un + Up)/2, agrees very well with 



predictions from many-body variational 



55l | and recent Dirac-Brueckner-Hartree-Fock 



(DBHF) [56| calculations. The underlying EOS has been also successfully tested against 
nuclear collective flow data in relativistic heavy-ion reactions at densities up to 5po (Po ~ 



0.16/m-3) |40 



571, 



58, 



591] . Also the strength of the symmetry potential estimated from the 
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> 



> 




FIG. 1: (Color online) Equation of state of stellar matter in /3-equilibrium. The upper panel shows 
the total energy density and lower panel the pressure as function of the baryon number density (in 
units of po). 

single-nucleon potentials via Ugym = {Un — Up) / {2a) at po agrees very well with the Lane 
potential extracted from nucleon-nucleus scatterings and (p, n) charge exchange reactions up 
to 100 MeV. The EOS outlined here has been applied recently to constrain the neutron-star 
radius [52|] with a suggested range compatible with the best estimates from observations. 

We show the EOSs used in this work in Fig. 1. The upper panel displays the total energy 
density, e, as a function of baryon number density and the lower frame shows total pressure. 
Predictions from Akmal et al. |60 | with the A18 + 5v + UIX* interaction (APR) and recent 



DBHF calculations (Bonn B) [33|,l34] with the Bonn B One-Boson- Exchange (OBE) potential 



are also shown. In addition to the pure nucleonic EOSs, we also include one hyperonic (Hyp) 
and one hybrid (Hyb) EOSs by the Catania group jsO], see e.g. Ref. 82]. The hyperonic EOS 



shown in Fig. 1 is an updated version of the one reported in Ref. 8l| • It has been calculated 



within the Brueckner-Hartree-Fock (BHF) approach using the Argone fig nucleon-nucleon 
(NN) potential, modified by nucleon three-body-forces according to the Urbana model, and 
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FIG. 2: (Color online) Neutron star mass, proton fraction, Yp, and symmetry energy, esym- The 
upper frame displays the neutron star mass as a function of baryon number density. The middle 
frame shows the proton fraction and the lower frame the nuclear symmetry energy as a function 
of density. (Symmetry energy is shown for the nucleonic EOSs only.) The proton fraction curve of 
the Hyb EOS is terminated at the beginning of the quark phase. The termination point is denoted 
by a "cross" character. 

with nucleon-hyperon interaction included. The hybrid EOS contains a BHF hadronic phase 
(n,p, e, yU, A, S) followed by a mixed phase {n,p,e, fi, A,T,,u,d, s), and a pure quark phase 
{u, d, s, e), calculated within the MIT bag model (with bag constant B = 90MeV fm"^ and 
150 MeV mass of the s-quark) [80]. The hadron-quark phase transition has been obtained 



by performing a Gibbs construction 



80|. Since, as demonstrated in Refs. [46|, [5^, only 



equations of state with x between -1 and have symmetry energies consistent with the 
isospin diffusion data and measurements of the skin thickness of ^^^Rb, we thus consider 



only these two limiting cases. Below the density of approximate 
of state shown in Fig. 1 are supplemented with a crust EOS 
for the low density regime. 



61 



y 0.07/m ^ the equations 



62l | which is more suitable 
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FIG. 3: (Color online) Mass-radius relation. The shaded region corresponds to the mass constraint 
from Hotan et al. 

TABLE 2: Maximum-mass neutron star models. The first column identifies the EOS. Remain- 
ing columns exhibit the following quantities for static configurations with maximum mass: total 
gravitational mass; radius; central baryon density. 



EOS 


Mmax{MQ) 


R{km) 




MDI(x=0) 


1.91 


9.89 


1.30 


MDI(x=-l) 


1.97 


10.85 


1.08 


Bonn B 


2.24 


10.88 


1.08 


APR 


2.19 


9.98 


1.14 


Hyp 


1.25 


9.96 


1.49 


Hyb (mixed phase) 


1.51 


9.14 


1.50 



Fig. 2 displays the neutron star mass (upper panel), the proton fraction (middle panel) 
and the nuclear symmetry energy (lower panel). The shaded region in the upper frame 
corresponds to the mass constraint by Hotan et al. [79]. From the neutron star models 
satisfying this constraint we observe that the proton fraction exceeds the direct URCA 
threshold for predictions from the x = —1 and Hyp EOSs. Here we recall that the fast 
URCA process proceeds only for Yp above 0.14 due to simultaneous conservation of energy 



and momentum 



631], causing rapid cooling of neutron stars. Hybrid stars (from the Hyb 



12 



EOS) cool rapidly through a quark direct URCA process (see Section V). 

The corresponding stellar sequences are shown in Fig. 3. Table 2 summarizes the gravi- 
tational masses, radii and central densities of the maximum-mass configurations. 



IV. CONSTRAINING THE CHANGING RATE OF THE GRAVITATIONAL 
CONSTANT G 

To constrain the hypothetical time variation of G, we consider stellar models constructed 
from different equations of state and calculate the neutron star surface temperature via the 
gravitochemical heating method. Since only EOSs with x between -1 and have symme- 



try energies consistent with terrestrial nuclear laboratory data [46|, [52], we consider these 
two limiting CcLSGS cLS representative of the possible range of neutron star structures. As 
demonstrated by Jofre et al. [1], if one assumes only slow Urea reactions in the neutron star 
interior, the stellar photon luminosity and surface temperature in the stationary state can 
be evaluated through Eqs. (13) and (15) respectively. In the most general case, however, 
when both direct and modified Urea reactions operate in the neutron star interior (as for 
the X = —1, Hyp and Hyb EOSs), the analytic expression (13) becomes a very poor ap- 
proximation [64]. Therefore, under this circumstance we calculate the surface temperature 
(and photon luminosity) in the quasi-equilibrium state by setting Eqs. (5), (8), and (9) to 
zero and solving them simultaneously. 

In Fig. 4 we show the neutron star stationary photon luminosity (upper panel) and steady 
surface temperature (lower panel) versus stellar mass, assuming \G/G\ = 4.5 x lO^^^yr^^. 
The value of G is chosen so that predictions from the x = EOS are just above the 90% 



confidence contour of Kargaltsev et al. 221], see Fig. 5. This upper limit is consistent with 



the one by Jofre at al. [if] under the same assumptions. Here we reiterate that we apply 



791 ] instead the one by van Straten et al. |23| used in 



the mass constraint by Hotan et al. 

Ref. 111]. (Both mass measurements partially overlap in the range (1.4 — 1.5)Mq.) We notice 
that predictions from the x = 0, APR and Bonn B EOSs all lie just above the observational 
constraints, with those from the x = and APR EOSs b eing very similar to each other. 



This observation has been already made in a previous work 52] in conjunction with a study 



of the neutron star radius and was interpreted in terms of the good agreement between the 
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corresponding symmetry energies up to about 5po (see also Fig. 2, lower frame). 

Figs. (6) and (7) display predictions assuming \G/G\ = 2.1 x lO^^^yr^^. In this case 
the value of G is chosen so that predictions from the x = — 1 EOS are just above the 
observational constraints. The surface temperature calculated with the x = —1 EOS is 
noticeably lower than predictions from EOSs allowing only modified Urea processes. This 
is due to opening of the fast Urea channel (see, Fig. 2, lower panel). When the neutron star 
mass becomes large enough for the central density to exceed the direct Urea threshold, tjhe 
surface temperature drops because of the faster relaxation toward chemical equilibrium 
The upper limit we provide here is more restrictive than the one reported by Jofre et al. 
in the case of fast cooling. 

Before we close the discussion in this section we also compare our results with recent 



Ti 



predictions from big bang nucleosynthesis (BBN) 70| ■ The stationary surface temperature in 
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FIG. 4: (Color online) Stationary photon luminosity (upper frame) and neutron star surface tem- 

G/G\ = A.5xW-^^yr~\ The shaded 



79[. The open (black) square charac- 



perature (lower frame) as functions of stellar mass, assuming 
region corresponds to the mass constraint form Hotan et al. 
ter denotes the maximum possible mass attained by the neutron star models constructed from the 
hyperonic EOS (Hyp). (The x = —1, Hyp and Hyb EOSs allow for fast cooling in the considered 
mass range. For the Hyb EOS, stellar cooling proceeds via a quark direct URCA process, see 
Section V.) 
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FIG. 5: (Color online) Neutron star stationary surface temperature for stellar models satisfying 



the mass constraint by Hotan et al. [79|. The solid lines are the predictions versus the stellar radius 
for the considered neutron star sequences. Dashed lines correspond to the 68% and 90% confidence 
contours of the black-body fit of Kargaltsev et al 
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IS 



chosen so that predictions from the x = EOS are just above the observational constraints. 
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that predictions from the x=-l EOS are just above the observational constraints. (The x = —1 
EOS allows direct Urea reactions in the considered mass range.) 

Fig. 8 is calculated assuming \G/G\ = 4x 10^^^yr~^ after Copi et al. 7o|. We observe that in 
this case the resulting surface temperatures are below the observational contour of Kargaltsev 
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FIG. 8: (Color online) Same as Fig. 5 but assuming \G/G\ = 4 x 10 ^^yr ^ after Copi et al. [70(]. 



et al. 221]. There are several possible reasons which could account for this discrepancy: (i) 



The time scales of both methods are quite different. While the approach applied by Copi et 
al. [3] probes the constancy of G over the age of Universe (see Table 1), the gravitochemical 



heating method is sensitive to variations in G over shorter periods 



69| . and falls closest 



to the second group of methods listed in Table 1. (ii) We do not consider the rotochemical 



heating [2^ in this paper. On the other hand, the rotochemical heating, which is believed to 
play an important role in the thermal evolution of millisecond pulsars, cannot explain alone 
the observations of Kargaltsev et al. [3] • Therefore, additional heating mechanisms need to 
be considered 2^. In fact, if taken together the roto- and gravitochemical heatings could, 
in principle, explain the empirical surface t emp erature of PSR J0437-4715, where a smaller 
upper limit for \G/G\ should be assumed [69|]. Under this circumstance, our predictions 
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and those of Copi et al. [TOj wou. 
account. As pointed out in Refs. 



d be relatively close, (iii) Superfluidity is not taken into 



69| , this would raise the surface temperature due to 



lengthening the time-scales needed to reach equilibrium. 



V. EFFECTS OF HYPERONS AND QUARKS 

Although the major goal of our work is to constrain the changing rate of G by applying 
an EOS with constrained symmetry energy, our analysis would be incomplete without a 
discussion of the possible effects of "exotic" states of matter in neutron stars. Namely, in 
what follows we discuss the effects of hyperons and quarks. 

Hyperons.- While at normal nuclear matter densities, p ~ Pq, neutron star matter con- 
sists of only nucleons and light leptons (e',/i~), at higher densities several other species of 
particles are expected to appear due to the rapid rise of the baryon chemical potentials with 
density [sill. The first particles to appear after the muons are the S~ and A° hyperons. 



Strange baryons appear mainly in reactions such as 51 1 



n + n <^==^ S +p (19) 
n + n <^ n + A° (20) 
n + A° +p (21) 

The equilibrium conditions for these reactions read 

2/i„ = /is- + yUp (22) 
l^n = I^AO (23) 

On the other hand, once the hyperons are present the following two chemical imbalances 
need to be introduced in addition to rjnpe and rinp/j, j3 

V2n^p = 2/i„ - /is- - /ip, (24) 
VnA = IJ'n- /^AO (25) 

Since processes (19) and (20) do not conserve strangeness, they can proceed only via weak 
interactions, while reactions (21) proceed via strong interactions. The above reactions have 



timescales at least several orders of magnitude shorter than those of beta processes [8J] and, 
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therefore, ?72nSp and ?7„a remain relatively small compared to T]npe and ?7„p^ J24]. Conse- 
quently, reactions (19)-(21) contribute negligibly to the total heat generation |24| . 

The importance of including Urea processes involving hyperons in the gravitochemical 
heating formalism, in addition to the nucleonic processes, can be assessed, for instance, by 



considering the following direct S Urea reactions 24 1 



^ n + l + u, (26) 
n + l S" + z/ (27) 

The net effect of these processes on Eqs. (8) and (9) is to enhance the lepton direct Urea 
rates, reducing the chemical imbalances {rjnpe, Vnpfi) and, in turn, the surface temperature. 
The overall correction is however small, because the direct Urea reactions with hyperons 
are at least a factor 5 weaker than the nucleon ones (see e.g. Prakash et al. jssi). For the 



case when only the modified Urea channel operates, Fernandez et al. [2J] showed that the 
correction to the surface temperature due to including several reactions involving hyperons 
is on the order of [Ln/{Ln + Lh)Y^'^^ ■, where Lat and Lh are the nucleon and hyperon Urea 
luminosities respectively. (This conclusion follows directly from Eq. (13).) Consequently, 



corrections due to hyperons in either direct or modified Urea processes can be neglected |24|. 

Yet, the inclusion of hyperons has a non- negligible effect on the EOS of neutron-star 
matter: Due to the appearance of additional degrees of freedom, the energy per particle 
is gre atly reduced and the resulting EOS is much softer than its purely nucleonic counter- 



part 8l|, l82] (see also Fig. 1). This effect is entirely due to the inclusion of hyperons as 



additional degrees of freedom and is observed in all current calculations independent of the 



adopted many-body approach and/or interaction (see e.g. [8^). The consequences for the 
neutron star global properties, namely radii and masses, are reflected in an overall very large 
mass reduction with respect to models of rape/x-stars (Fig. 3). 

Thus, in light of the above considerations, to investigate the impact of hyperons on the 
changing rate of G we apply a hyperonic EOS (see Section III) and the gravitochemical 
heating formalism without including the imbalances (24) and (25) (i.e., we do not consider 
openings of the hyperonic channels for additional heat emission). Rather, we aim to identify 
the effect of the global neutron star properties (masses and radii) from hyperonic EOS on the 
time variations of G in a qualitative way. In calculating the surface temperature we adopt a 
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FIG. 9: (Color online) Same as Fig. 5 but assuming \G/G\ = 3.6 x 10 ^^yr ^. The value of G is 
chosen so that predictions from the hyperonic EOS (Hyp) are just above the empirical constraints. 

fast cooling scenario since including of hyperonic degrees of freedom increases proton fraction 
above the threshold for fast nucleonic Urea processes for neutron stars in the considered mass 
range (see Fig. 2 upper and middle frames). We show our predictions for \G/G\ based on 
the hyperonic EOS (Hyp) in Fig. 9. The value \G/G\ = 3.6 x lO^^^yr^^ reassures that 
results from the Hyp EOS are (just) above the empirical contours and is comparable with 
the one predicted with the x = —1 EOS (which also allows for fast Urea reactions). Here 
a few comments are in place: First, we should mention that due to the "softness" of the 
Hyp EOS, the maximum mass attained by the neutron star models in this case is much 
smaller than the corresponding masses predicted by purely nucleonic EOSs (see Fig. 3 and 
Table 1). The maximum mass is about 1.25Mq and falls closer to the lower limit of the 
allowed mass range for PSR J0437-4715. In fact, if one adopts the constraint by van Straten 



et al. |23| {MpsR = 1.58 ± O.ISM©), Mmax would be well below the allowed range. Also, as 



pointed out by Burgio et al. 82], a maximum mass less than I.3M0 is in contradiction with 
the most precisely measured value of the Hulse- Taylor pulsar mass, PSR 1913+16, which 
amounts to 1.44Mq; Second, hyperonic stars are much more compact than rape/i-stars. With 
the Hyp EOS the maximum mass is achieved at pc ~ 1.48/m~^ (see Table 2). The rapid 
rise of the central density with increasing/decreasing neutron star mass/radius causes a 
steep growth in the predicted surface temperature (see e.g.. Fig. 9); Finally, in the context 
of the gravitochemical heating of interest here, the only appreciable effect of hyperons is 
that their presence could facilitate the fulfilment of the direct Urea conditions, as in the 
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case of reactions involving A° [83]. Additionally, the rise of the proton fraction due to the 
appearance of E~ hyperons can shift the threshold of the nucleon direct Urea process to 
lower densities |85|. These considerations, together with the Hyp EOS predictions in Fig. 9, 



support the findings of Reisenegger 86], namely that the EOSs do not give very different 
results for the surface temperature if one only compares among those which either do or do 
not allow for direct Urea reactions. 

Quarks.- The existence of quark matter in the interior of neutron stars is one of the 
major issues in the physics of these compact astrophysical objects. Here we discuss how the 
appearance of deconfined quarks would affect a possible time variation of G. For our analysis 
we choose a hybrid EOS (Hyb) derived by the Catania group (for details, see Section III 
and the references therein). 

Quark matter is assumed to consist of deconfined u, ci, and s quarks, and a small fraction 
of electrons [87]. At energies relevant to neutron stars, the u and d quarks are treated 
as massless particles, while the s quark is moderately relativistic and has a non-negligible 
mass jsi]- The most important processes among the constituents of the quark-matter system 



are the direct Urea reactions with u and d quarks [Ml] 



d u + e + Ue, (28) 
u + e d + Ue (29) 

In addition, the following direct Urea reactions involving the s quark are possible: 

s — > M + e + z/g, (30) 
u + e —>■ s + Ve (31) 

These /3-processes, however, do not conserve strangeness and, generally, yield much smaller 



emissivity relative to reactions (28) and (29) 51|- The chemical equilibrium conditions imply 



IJ'd = IJ-s = f^u + fJ'e (32) 

A rigorous treatment of quark matter in the gravitochemical formalism would require 
major adjustments/extensions of the framework to incorporate, for instance, the imbalances 
rjued and rjues, and the corresponding quark emissivities. Since this task is beyond the main 
focus of this work, we take a much simpler approach which should, at minimum, hint on 
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FIG. 11: (Color online) Same as Fig. 5 but assuming \G/G\ = 1.6 x 10 ^^yr ^. The value of G is 
chosen so that predictions from the hybrid EOS (Hyb) are just above the empirical constraints. 

the effect of quarks relative to the npe/x-matter. Similarly to what is done in the case of 
hyperonic stars, we provide a (qualitative) estimate based on the mass-radius correlation of 
hybrid stars, rather than their detailed composition. Because the triangle inequality 

fc^ < fc]^ + k'p, I = {d, s} (33) 
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among the Fermi momenta of species u, d, s and e, is fulfilled for reactions (30) and (31) 
but is not for (28) and (29) (see Fig. 10), cooling of hybrid stars (from the Hyb EOS) can 
proceed through quark direct Urea processes involving the s quark. Therefore we assume 
a fast-cooling scenario. Our predictions for \G/G\ based on the Hyb EOS are shown in 
Fig. 11. We notice that the value \G/G\ = 1.6 x 10~^^yr~^ (chosen so that predictions 
from the Hyb EOS are above the empirical data) falls closest to the one derived by Jofre 
et al. [l| under the modified Urea assumption (and also our predictions from the x = —1 
and Hyp EOSs). Since the yield of reactions (30) and (31) is less certain [5l| and because 
we do not treat the inclusion of quarks explicitly, we should keep in mind that the present 
result is more qualitative than quantitative and therefore more refined calculations might 
be necessary to draw more definite conclusions. On the other hand, we observe a steep 
rise of surface temperature with stellar mass (red-shifted radius and/or central density), see 
e.g. Fig. 11. The rapid temperature increase is attributed to the fast growth of central 
density with mass (Fig. 10, upper panel), and has been already observed in conjunction 
with the analysis of hyperonic stars (see above). Finally, we notice that hybrid stars have, 
generally, smaller radii than those of npe/i-stars (Fig. 3). In fact, the radii predicted with 
the Hyb EOS employed here are below the latest neutron star radius constraints j52] (see 
also Fig. 3). Yet, models yielding quite stiff quark matter EOSs, and, in turn, hybrid/quark 
star radii (masses) consistent with the latest observations [66\, do exist [88!|. Given the 
large present uncertainties in modeling the quark phase of matter, our results provide a 
reasonable qualitative assessment of the changing rate of G, based on the gravitochemical 
heating method and a hybrid EOS. 



VI. SUMMARY 



In summary, to test the Dirac hypothesis and constrain the changing rate of the grav- 
itational constant G due to the continuous expansion of the Universe, we have calculated 
the neutron star surface temperature through the gravitochemical formalism introduced by 
Jofre et al. |l| applying several EOSs spanning from pure nucleonic matter to quark matter. 
One of the nucleonic EOSs has a density dependence of the symmetry energy constrained by 
terrestrial nuclear laboratory data that just became available recently. Using the "softer" 
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symmetry energy {x = 0) consistent with both the isospin diffusion and the "^^^Ph neutron 
skin data, we obtain an upper hmit \G/G\ < 4.5 x lO^^^yr^^ compatible with that obtained 
by Jofre et ah l| under the same assumptions. This is mainly because the density depen- 
dence of the symmetry energy with x = turns out to be very close to that with the APR 
EOS they also used. The "stiffer" symmetry energy (x = —1) EOS yields an upper limit 
\G/G\ < 2.1 X 10~^^yr~^, an order of magnitude lower than the one derived by Jofre et al. |l] 
when both direct and modified Urea reactions operate. Since both the x = and x = — 1 
equations of state have symmetry energies consistent with the empirical nuclear data, we 
cannot rule out the results from the direct Urea scenario. In the case of fast cooling, our 
result provides even tighter upper limit than the one reported in Ref. Ij. 

The gravitochemical heating mechanism has the potential to become a powerful tool for 
constraining gravitational physics. Since the method relies on the detailed neutron star 
structure, which, in turn, is determined by the EOS of stellar matter, further progress in 



our understanding of properties of dense, neutron-rich matter wil 
effective. Precise astrophysical observations such as those by Ozel 



make this approach more 



66| , Hessels et al. [65| , and 



Kaaret at al. 67|] together with future heavy-ion experiments with high energy radioactive 



beams 



68l | will allow us to set more stringent constraints on the EOS of dense neutron-rich 



matter, and on the possible time variation of the gravitational constant G. 
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